使用C语言和Python基于MPI并行计算PI的值

#使用C语言和Python基于MPI并行计算PI的值 | 来源: 网络整理| 查看: 265

浏览量： 1,303 (在苹果系统下，如果文章中的图片不能正常显示，请升级Safari浏览器到最新版本，或者使用Chrome、Firefox浏览器打开。)

圆周率PI是一个很神奇的数字，自古以来，包括数学家在内的很多人都曾使用过各种各样的算法去估算PI的真实值，并且都取得了一定的成就。古巴比伦人使用3.125作为PI的近似值，约公元前1700年的古埃及人则提出PI=3.1604，中国的祖冲之(430-501)则使用355/113作为近似值，使得PI值精确到了7位数。随着计算机的问世，以及科学技术发展的需要，PI的近似值目前精确位数早已突破万亿位。PI值除了有其每一位、每两位、每三位都符合均匀分布的统计规律特性之外，还可以用来检测计算机硬件的可靠性，而且，也可以用来入门并行计算。

关于PI的计算的各种各样的方法，以及其发展历程，可以参考知乎：

用计算机算圆周率，是个怎样的过程？

在这里，我将使用两种方法来实现，并且分别再使用C语言和python实现一遍。

方法1. 无穷级数法 J. Gregory(1638-1675)和G. W. Leibniz(1646-1716):

PI / 4 = 1 – 1/3 + 1/5 – 1/7 + ……

C语言

//mpic++ 1.cpp -o 1.out //mpiexec -n 32 ./1.out #include #include #include "mpi.h" double calculate_part(long long start, long long step); int main(int argc, char *argv[]) { int my_rank, comm_size; clock_t t_start, t_finish; double t_duration; long long N = (long long)1024 * (long long)1024 * (long long)64 * (long long)1; MPI_Init(&argc, &argv); MPI_Comm_size(MPI_COMM_WORLD, &comm_size); MPI_Comm_rank(MPI_COMM_WORLD, &my_rank); long long step = (long long)(N/comm_size); long long start = my_rank * step; t_start = clock(); double value = calculate_part(start, step); double result = 0.0; if (my_rank == 0) { result += value; for(int i = 1; i < comm_size; i++) { MPI_Recv(&value, 1, MPI_DOUBLE, i, 0, MPI_COMM_WORLD, MPI_STATUS_IGNORE); result += value; } t_finish = clock(); t_duration = (double)(t_finish - t_start) / CLOCKS_PER_SEC; printf("pi is : %0.15f\n",result * 4); printf( "%f seconds\n", t_duration ); } else { MPI_Send(&value, 1, MPI_DOUBLE, 0, 0, MPI_COMM_WORLD); } MPI_Finalize(); return 0; } double calculate_part(long long start, long long step) { double sum = 0.0; int flag = 1; for(long long i = 0; i < step; i++) { if(start%2!=0) flag = -1; else flag = 1; sum += flag * (1/(double)(2*start+1)); start ++; } return sum; }

运行测试

启动进程节点数计算任务量 PI值结果时间开销 (s) 1 1024*1024*64 3.141592653471873 28.302326 2 1024*1024*64 3.141592653471685 14.391914 4 1024*1024*64 3.141592653471855 8.191455 8 1024*1024*64 3.141592653471702 4.418357 16 1024*1024*64 3.141592653472924 2.392712 32 1024*1024*64 3.141592653472887 1.805563

加速比(1,2) = 1.9665 加速比(2,4) = 1.7569 加速比(4,8) = 1.8540 加速比(8,16) = 1.8466 加速比(16,32) = 1.3252

python

from mpi4py import MPI import time comm = MPI.COMM_WORLD rank = comm.Get_rank() size = comm.Get_size() def calculate_part(start, step): sum=0.0 flag=1 for i in range(0,step): if(start % 2 != 0): flag=-1 else: flag=1 sum+=flag * (1/(2*start+1)) start +=1 return sum N = 1024*1024*64 step = N // size start = rank * step t0=time.time() value = calculate_part(start, step) result = 0.0 if rank == 0: result += value for i in range(1,size): value = comm.recv(source=i, tag=0) result += value t1=time.time() print('PI is : ',result * 4) print('time cost is', t1 - t0, 's') else: comm.send(value, dest=0, tag=0) 方法2. 蒙特卡洛模拟

C语言

//mpic++ -std=c++11 2.cpp -o 2.out //mpiexec -n 32 ./2.out #pragma GCC diagnostic error "-std=c++11" #include #include #include #include #include #include "mpi.h" using namespace std; double calculate_part(int step); int main(int argc, char *argv[]) { int my_rank, comm_size; clock_t t_start, t_finish; double t_duration; int N = 1024 * 1024 * 64; MPI_Init(&argc, &argv); MPI_Comm_size(MPI_COMM_WORLD, &comm_size); MPI_Comm_rank(MPI_COMM_WORLD, &my_rank); int step = (int)(N/comm_size); t_start = clock(); double value = calculate_part(step); double result = 0.0; int count = 0; if (my_rank == 0) { count += value; for(int i = 1; i < comm_size; i++) { MPI_Recv(&value, 1, MPI_DOUBLE, i, 0, MPI_COMM_WORLD, MPI_STATUS_IGNORE); count += value; } result = 4.0 * ((double)count / (double)N); t_finish = clock(); t_duration = (double)(t_finish - t_start) / CLOCKS_PER_SEC; printf("pi is : %0.15f\n",result); printf( "%f seconds\n", t_duration ); } else { MPI_Send(&value, 1, MPI_DOUBLE, 0, 0, MPI_COMM_WORLD); } MPI_Finalize(); return 0; } double calculate_part(int step) { std::default_random_engine random(time(NULL)); std::uniform_real_distribution dis(0.0, 1.0); int count = 0; for(int i = 0; i < step; i++) { double x = dis(random); double y = dis(random); //int r = pow(x*x+y*y,0.5); double r = x*x+y*y; if(r AI柠檬博主正在阿里云上出售域名 “y403.com”，感兴趣就快去看看吧打赏赞(2)微海报分享

【本文地址】

公司简介

联系我们